#include    <RootFinder.h>

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const UInt  RootFinder::MAX_ITERATIONS = 20;

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void RootFinder::FindRoots
    (
    Function&       func,
    const Real&     from,
    const Real&     to,
    RealList&       roots
    )
{
//  Initialize the return result set
    roots.clear();

//  Get singularity points for the function
    RealList        singularities;
    func.GetSingularities(from, to, singularities);

//  Add the "from" point into the head of the list and "to" point to its tail
    singularities.push_front(from);
    singularities.push_back(to);

//  Now, loop through the list, interval-by-interval, and search for a root
//  at each interval

    RealListCI      intervalFrom    = singularities.begin();
    RealListCI      intervalTo      = intervalFrom;
    ++intervalTo;

    while (intervalTo != singularities.end()) {
        Real    xLeft   = *intervalFrom;
        Real    xRight  = *intervalTo;

        ++intervalFrom;
        ++intervalTo;

//      Adjust the interval to step away from possible asymptotes
        if (RootFinder::AdjustInterval(func, xLeft, xRight) == AOK) {
            Real    root;
            if (RootFinder::FindRoot(func, xLeft, xRight, root) == AOK)
                roots.push_back(root);
        }
    }
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Status RootFinder::FindRoot
    (
    Function&       func,
    const Real&     from,
    const Real&     to,
    Real&           root
    )
{
    Real    xLeft   = from;
    Real    yLeft   = 0;
    func.GetValue(xLeft, yLeft);
    Int     signLeft = RootFinder::Sign(yLeft);

    Real    xRight  = to;
    Real    yRight  = 0;
    func.GetValue(xRight, yRight);
    Int     signRight = RootFinder::Sign(yRight);

//  A very special case when both boundaries are roots and yet
//  they have different values. This should never happen but
//  for genericness we have to support that.
    if ((signLeft == signRight) && (signLeft == 0) && !(xLeft == xRight)) {
        return ERROR;
    }

//  Special case when left boundary is a root
    if (signLeft == 0) {
        root = xLeft;
        return AOK;
    }

//  Special case when right boundary is a root
    if (signRight == 0) {
        root = xRight;
        return AOK;
    }

//  See if root can be found. 
    if (signLeft == signRight) {
        return ERROR;
    }

    Real        xMiddle     = 0;
    Real        yMiddle     = 0;
    Int         signMiddle  = 0;
    for (UInt i = 0; i < MAX_ITERATIONS; ++i) {
        xMiddle = (xLeft + xRight) / Real_TWO;
        func.GetValue(xMiddle, yMiddle);
        signMiddle = RootFinder::Sign(yMiddle);
        if (signMiddle == signLeft) {
            xLeft = xMiddle;
            continue;
        }

        if (signMiddle == signRight) {
            xRight = xMiddle;
            continue;
        }

        xLeft  = xMiddle;
        xRight = xMiddle;
        break;
    }

    root = (xLeft + xRight) / Real_TWO;

    return AOK;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Int RootFinder::Sign(const Real& x)
{
    if (x > Real_ZERO)
        return 1;

    if (x < Real_ZERO)
        return -1;

    return 0;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Status RootFinder::AdjustInterval
    (
    Function&       func,
    Real&           from,
    Real&           to
    )
{
//  Let's set the initial offset at 1/3 of the original interval length
//  It has to be less than 1/2 of it, though.
    Real    offset = (to - from) / Real_THREE;

    for (UInt i = 0; i < MAX_ITERATIONS; ++i) {
        Real    xLeft   = from + offset;
        Real    yLeft   = 0;
        func.GetValue(xLeft, yLeft);
        Int     signLeft = RootFinder::Sign(yLeft);

        Real    xRight  = to - offset;
        Real    yRight  = 0;
        func.GetValue(xRight, yRight);
        Int     signRight = RootFinder::Sign(yRight);

//      Special case when left boundary is a root
        if (signLeft == 0) {
            offset /= Real_TWO;
            from    = xLeft - offset;
            to      = xLeft + offset;
            return AOK;
        }

//      Special case when right boundary is a root
        if (signRight == 0) {
            offset /= Real_TWO;
            from    = xRight - offset;
            to      = xRight + offset;
            return AOK;
        }

//      See if root can be found on the current interval
        if (signLeft != signRight) {
//          Yes, it can!
            from    = xLeft;
            to      = xRight;
            return AOK;
        }

//      Let's make the offset smaller...
        offset /= Real_TWO;
    }

//  If we reached that far, we could not stretch the interval
//  big enough to enclose a root.
    return ERROR;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// eof
